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ABSTRACT 


This  thesis  investigates  and  demonstrates  the  workability  of  the  time-reversed 
process  for  radar  imaging  applications,  particularly,  for  bi- static  or  multi- static  radars. 
One  benefit  of  the  time-reversed  process  is  its  ability  to  reduce  the  calculation  to 
determine  the  targets’  shape.  The  finite- difference- time- domain  (FDTD)  method  is  used 
to  demonstrate  the  time-reversed  process. 

Following  an  overview  and  description  of  the  principles  of  the  time-reversed 
process,  the  FDTD  method  is  applied  to  the  wave  equation  and  the  time  reversed-process 
in  2-D  space.  The  FDTD  numerical  model  is  developed  and  used  for  producing 
fundamental  examples  on  conducting  targets.  The  examples  reveal  that  the  time-reversed 
process  can  be  employed  for  radar  imaging  within  certain  constraints.  Finally, 
conclusions  regarding  the  time-reversed-process  are  presented  and  recommendations  for 
future  research  are  provided. 
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1.  INTRODUCTION 


A.  OVERVIEW 

Bi- static  or  multi- static  radar  has  been  recognized  as  an  effective  tool  for 
detecting  and  tracking  low  observable  targets,  such  as  stealth  aircraft  [Ref.l].  The 
development  of  these  systems  emphasized  measurement  of  the  target’s  precise  location  in 
range  and  angle.  However,  in  general,  they  are  not  capable  of  identifying  the  shape  of  a 
target.  In  many  situations,  this  information  may  be  important  to  make  an  appropriate 
decision  or  to  distinguish  a  friendly  target  from  the  enemy.  Because  of  its  importance, 
radar  imaging  has  been  a  topic  of  intense  research  for  the  past  several  decades,  being 
driven  by  the  need  to  identify  friend  from  foe  in  future  radar  systems. 

At  present  some  useful  radar  techniques  are  available,  such  as  SAR/ISAR,  which 
produce  very  detailed  information  of  a  target.  These  techniques  require,  however,  a 
relatively  long  process  in  obtaining  the  information.  Therefore,  more  timely  processing 
techniques  are  desired  in  many  situations. 

To  obtain  information  relating  target  identity  to  the  received  radar  signature, 
inverse  scattering  problems  have  to  be  solved.  Solutions  to  the  forward  scattering 
problem,  which  predict  scattering  from  known  objects,  are  diverse  and  well  known. 
Once  the  scattered  wave  is  known  imphcitly  or  exphcitiy,  determining  a  target’s  shape  is 
possible,  in  theory,  by  inverting  the  forward  scattered  solution  process.  Many  methods 
for  inverse  scattering  have  been  investigated.  For  example,  Colton  and  Kress  [Ref.2] 
introduced  inverse  scattering  theory  in  which  Colton,  Giebermann  and  Monk  [Ref.3] 
provided  numerical  examples  of  obtaining  the  shape  of  an  obstacle  in  three  dimensions 
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using  time-harmonic  incident  and  scattered  acoustic  waves.  These  methods  were  induced 
by  solving  the  frequency-  domain  Helmholtz  equation. 

Some  useful  methods  for  solving  scattering  problems  in  the  time-  domain  exist, 
such  as  time- domain  physical  optics,  time- domain  integral  equations  and  finite- difference 
time- domain  (FDTD).  These  methods  can  be  used,  in  theory,  to  determine  or  estimate 
the  shape  of  a  scattering  obstacle.  Although  there  are  many  pubhcations  available  on  the 
inverse  scattering  time-domain  methods,  these  deal  mainly  with  an  electric  or  magnetic 
field  integral  equation  [Ref. 4].  The  integral  equation  approach  requires  enormous 
computational  resources  to  determine  the  shape  of  a  complex  target  and  requires  a  very 
high  signal  to  noise  ratio  in  the  received  scattered  signatures. 

This  thesis  investigates  a  new  technique  for  estimating  the  location  and  shape  of 
one  or  more  targets  by  processing  scattering  signature  information  captured  by  sensors  on 
the  perimeter  of  a  region.  A  reversed- time  solution  within  the  region  is  performed  using 
the  FDTD  method.  Such  a  method  has  the  potential  for  employment  by  bi- static  or 
multi- static  radar  systems. 

Fink  [Ref.5]  introduced  the  concept  of  “Time  Reversed  Acoustic  Imaging” 

associated  with  wave  field  propagation.  Fink’s  colleagues  [Ref.  6,  7,  8,  9]  have  worked 

with  this  concept.  The  basic  idea  of  this  process  is  based  on  an  elementary  fact  of  time- 

reversal  invariance  of  the  wave  equation  in  a  lossless  medium.  Briefly,  the  process  occurs 

when  the  wave  field  from  an  obstacle  or  a  target  is  propagated  and  is  captured  on  the 

surrounding  surface  using  a  number  of  transducers.  These  transducers  combine  the 

functions  of  microphone  and  loudspeaker,  emitting  a  time-reversed  replica  of  the 

received  signature  from  each  transducer.  These  emissions  generate  a  time-reversed  field 
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within  the  enclosed  region,  which  coUapses  upon  the  target.  In  other  words,  the  wave 
field  is  recorded  on  the  boundary  and  is  reradiated  as  an  appropriate  reversed  field.  The 
reversed  field,  then,  focuses  on  the  original  source.  This  process  may  be  extended  to 
electromagnetic  waves. 

The  ‘Time  Reversed  Imaging”  technique,  which  is  the  subject  of  study  in  this 
thesis,  uses  the  “Time  Reversed  Acoustics”  concept.  This  concept  focuses  on  wave  field 
propagation,  in  particular  wavefronts.  This  imphes  that  the  field  is  not  always  present  at 
the  target  location,  which  may  reduce  computational  time.  Moreover,  focusing  on 
wavefronts  enables  this  concept  to  be  used  for  a  wide  frequency  range  and  for  weak 
scattering  cases.  In  general  terms,  this  concept  seems  simple  and  logical,  but  historically 
it  has  not  been  used  for  electromagnetic  waves. 

B.  THESIS  OBJECTIVE 

This  thesis  investigates  apphcation  of  the  time-reversed  process  to 
electromagnetic  waves  and  employment  for  bi- static  or  multi- static  pulsed  radar 
apphcations.  The  thesis  develops  needed  numerical  modehng  and  simulation  using  the 
FDTD  method,  which  is  an  efficient  way  to  represent  this  process  for  vector  fields. 

In  Chapter  2,  the  general  principle  of  the  time  reversal  process  is  introduced. 
Chapter  3  focuses  on  the  numerical  representation  and  description  of  modeling  and 
simulation.  Chapter  4  presents  numerical  examples  using  MATLAB.  Finally,  Chapter  5 
summarizes  and  concludes  this  study  with  suggestions  for  further  research. 
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n.  GENERAL  PRINCIPLE  OE  TIME-  REVERSED  PROCESS 


This  chapter  presents  an  overview  of  the  time -reversed  process.  Details  can  be 
found  in  Reference  8. 


A.  TIME  REVERSAL  OF  ELECTRIC  FIELDS 

Maxwell’s  equations  induce  the  wave  equation  for  the  electric  field  in  lossless 

homogeneous  media.  In  free  space,  this  is  written  as 


V"E(r,t) 


1  a"E(r,t) 
Z  at" 


(1) 


where  V  represents  the  Laplacian  operator  with  respect  to  the  spatial  r-  coordinates  and 
c  represents  the  phase  velocity  in  free  space. 


Looking  at  equation  (1),  it  contains  only  a  second  order  time-derivative  operator. 
It  follows,  therefore,  that  if  ^(r,  t)  is  a  solution,  then  ^(r,  -t)  is  also  a  solution.  That  is,  the 
wave  equation  is  unchanged  under  a  time-reversal  transform  if  there  is  no  absorption 
during  propagation  in  the  medium.  This  property  is  the  starting  point  of  the  time- 
reversal  principle.  In  order  for  the  time-reversal  invariance  to  remain  valid,  no  loss  or 
absorption  during  propagation  is  assumed. 

The  special  property  of  time  reversal  has  been  observed  by  Stokes  [Ref.  9]  in  the 

classical  experiment  of  reflection  and  transmission  of  a  plane  wave  between  two  different 

media.  Considering  an  incident  plane  wave  of  amphtude  Eq\  propagating  from  medium  1 

to  medium  2,  we  can  observe  the  amphtude  of  the  reflected  field  (Eor)  and  the  transmitted 

field  (Eot)-  Here,  introducing  r  as  the  amphtude  reflection  coefficient  and  t  as  the 
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amplitude  transmission  coefficient  leads  to  ^or  =  r  Eqi  and  £'ot=  t  E’oi.  Following  this,  a 
solution  of  the  wave  equation,  ^(r,  t),  results  from  three  plane  waves.  Figure  1  indicates 
this  situation. 

In  the  case  of  no  absorption,  a  wave  must  be  reversible  in  accordance  with  the 
reciprocity  theorem  leading  to  the  principle  of  reversibihty.  Thus,  the  condition  of 
Figure  2  must  also  be  physically  possible.  Then,  the  time-  reversed  solution,  ^(r,  -t),  can 
be  described  by  Eoi,  Eq^  and  ^ot.  Accordingly,  when  examining  the  situation  in  Figure  3, 
there  are  two  incident  waves  of  amplitude:  xE^i  and  Coi-  One  wave  whose  amplitude  is 
ffioi  is  both  reflected  and  transmitted  at  the  interface.  Letting  r’  and  t’  be  the  amplitude 
reflection  and  transmission  coefficient,  respectively,  for  a  wave  incident  from  medium  2 
to  medium  1,  the  reflected  portion  is  tr’E'oi  and  the  transmitted  portion  is  tt’E'oi.  Similarly, 
the  incoming  wave  whose  amphtude  is  Eoi  sphts  into  both  the  amphtude  rr  Eoi  and  rt  Eo{. 
Using  the  superposition  principle,  it  follows  that 

tt’£oi+rr£oi  =  £'oi  (2) 

tr’E’oi-i- rt£'oi=0  (3) 

Therefore, 

tt’  =  1  -  r^  (4) 

r’  =  -r  (5) 

If  the  interface  between  medium  I  and  medium  2  is  a  perfect  conductor,  t  and  t’ 
are  zero.  This  results  in  either  1,  r’=  -1  or  r=  -1,  r’=l.  For  simphcity,  this  thesis 
considers  only  targets  or  obstacles  whose  surfaces  have  the  condition  of  a  perfect 
conductor.  Furthermore,  it  is  important  to  note  that  the  two  relationships  written  above 
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are  only  valid  if  the  reflected  and  transmitted  plane  waves  are  able  to  propagate  without 
attenuation,  which  knphes  that  they  have  a  real  wave  number.  The  incident  field  we 
consider  in  this  thesis  does  not  contain  any  evanescent  waves  that  cannot  be  time- 
reversed. 


Figure  1 .  Reflection  and  transmission  of  a  plane  wave  along  the  interface  between 
medium  1  and  medium  2 
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Figure  2.  Time-  reversal  of  Figure  1 


Figure  3.  Explanation  of  time-  reversal  for  reflection  and  transmission 
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B.  PRACTICAL  CONSIDERATION 

The  initial  conditions  of  a  target  or  an  obstacle  and  the  boundary  conditions 

determine  a  unique  solution,  ^(r,  t)  to  the  wave  equation  in  (1).  In  the  time-reversed 
operation,  the  solution  ^(r,  -t)  is  given  by  modifying  the  initial  condition  and  the 
boundary  condition.  In  practice,  however,  because  every  physical  phenomenon  requires 
causahty,  ^(r,  -t)  is  not  a  valid  solution.  Therefore,  we  consider  ^(r,  T-t)  under  the 
limited  time  interval  T,  where  T  is  sufficiently  advanced  in  time  so  that  ^(r,  t)  is  regarded 
as  zero  for  t  >  T.  To  assume  otherwise  requires  “initial  condition”  knowledge  of  the 
fields  ^(r,  t)  in  the  whole  three-dimensional  volume  daring  the  time  interval  T  in  order  to 
generate  the  time-reversed  solution  ^(r,  T-t).  A  more  realistic  way  to  generate  the 
time-reversed  solution  is  to  use  the  advantage  of  Huygens’  principle.  Based  on  this 
principle,  the  time-reversed  operation  in  the  three-dimensional  volume  requires  time- 
reversed  boundary  conditions  on  an  enclosing  two-dimensional  surface.  Using  this 
approach,  focusing  on  a  target  can  be  described  in  the  following  way. 

During  the  recording  step,  the  target  within  the  volume  surrounded  by  a  finite 
number  of  receivers  generates  a  field,  ^(r,  t),  which  produces  an  expanding  wavefront. 
The  receivers  sample  the  field  at  locations  on  the  enclosing  surface  and  have  the 
capability  to  measure  the  field  without  disturbing  the  propagation  of  the  field. 
Therefore,  the  field  is  propagated  in  a  free  unbounded  space.  During  the  time-reversed 
reconstruction,  the  target  behaves  as  a  passive  source  or  is  ignored.  Each  receiver  then 
generates  the  field,  ^(r,  T-t),  that  corresponds  exactly  to  the  time- reversal  of  the 
corresponding  field  measured  during  the  recording  step.  A  time-reversed  field  back- 
propagates  inside  the  recording  surface  and  is  focused  on  the  initial  target  position. 

Figures  4  and  5  illustrate  the  recording  and  reconstmction  steps,  respectively. 
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Figure  4.  Recording  step  of  the  time-  reversal  process  with  a  closed  surface 
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Figure  5.  Reconstruction  step  of  the  time -reversal  process  with  a  closed  surface 
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m.  NUMERICAL  MODELING  OE  TIME  -REVERSED  PROCESS 


This  chapter  presents  an  overview  of  the  finite- difference  time- domain  (FDTD) 
numerical  model  used  for  demonstrating  the  time-reversed  process  in  2-D  space.  The 
FDTD  method  is  an  appropriate  way  to  represent  this  process,  applying  approximations 
for  the  space  and  time  partial  derivatives  within  the  wave  equation.  Space  and  time 
discretizations  are  selected  to  ensure  the  numerical  stabihty  of  the  algorithm. 

For  the  region  defined,  the  wave  equation  is  solved  subject  to  initial  value  and 
boundary  value  conditions  according  to  electromagnetic  theory.  This  is  the  same 
procedure  apphed  to  the  FDTD  method.  Overall,  the  FDTD  is  a  time -marching  procedure 
discretely  simulating  the  continuous  fields  that  solve  the  wave  equation.  The  solution  is 
subject  to  the  initial  and  boundary  conditions.  At  each  time  step,  the  system  of  equations 
is  used  to  update  the  field  components  and  is  fuUy  explicit  so  that  setting  up  or  solving 
linear  equations  is  not  required.  Time  stepping  is  continued,  allowing  one  to  observe 
electromagnetic  wave  propagation  and  scattering.  Details  on  the  FDTD  can  be  found  in 
many  pubhcations.  Some  of  these  are  fisted  in  References  10  and  11. 

The  next  section  presents  numerical  modeling  of  the  wave  equation  followed  by 
the  forward  and  reverse  processes,  and  is  based  on  Reference  12.  For  simplicity,  a  3-D 
(2-dimentional  plus  time)  rectangular  region  with  E^(x,y,t)  (TMz  waves)  is  assumed 
throughout  this  thesis. 
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A.  THE  FDTD  SOLUTION  OF  THE  WAVE  EQUATION 

This  section  briefly  provides  details  on  how  the  wave  equation  can  be  solved  with 

the  FDTD  method. 


Consider  solving  the  wave  equation  in  a  rectangular  region,  0=  x 

=  a,  0=  y  =  b  and 

t  =  0,  as  depicted  in  Figure  6, 

V^E^  (X,  y,  0  =  0 

(6) 

with  initial  conditions 

£'^(x,y,0)  =  /(x,y) 

(7) 

3E 

"  (x,y,0)=  g(x,y) 
dt 

(8) 

and  boundary  conditions. 

E^(0,y,t)  =  ei(y,t),  E^{a,y,t)  =e^{y,t) 

E^(x,0,t)  =  e^(x,t),  E^(x,b,t)  =  e^(x,t) 

(9) 

Here,  the  assumption  is  that  a  =  Mh  and  b=  Nh,  where  ?v:=?y=h  is  an  equal  grid 
spacing  in  x  and  y  with  M  and  N  as  arbitrary  integers.  Letting  ?t  be  the  time  step  and  T 
be  the  desired  time  interval,  the  grid  positions  are  given  by 

=  (m-l)h,  y„  =  (n-l)h  and  =  (p-l)At  (10) 

where  m  =1,M,  n  =  1,V  and  p=\,T  . 

We  deflne  {m,n,  p)  to  be  the  discrete  space-time  sample  of  p) . 

Applying  the  finite  difference  approximations  in  space  and  time  at  the  point 
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Equation  (14)  implies  that  we  can  predict  at  the  point  yn,  tp+i)  from  six  previous 
points  in  the  grid  as  shown  in  Figure  7. 

By  introducing  2-D  spatial  arrays  Ep  =Ez(:  ,  :  ,  p),  where  Ep  is  an  MxN  array  at 
time  step,  p,  we  can  rearrange  equation  (14)  to  the  form, 
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where 
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are  respectively  MxM  and  NxN  unit  off-diagonal  arrays.  Equation  (14)  is  the  time- 
marching  equation  for  non-boundary  terms.  To  insure  the  stabihty  of  solutions,  the  next 
condition  must  satisfy  the  2-D  Courant  condition: 


cAt  ^  1 

h  ~  42 


(15) 


At  p  =  1,  the  initial  condition  (7)  is  used  and  the  value  at  p  =  0  is  established  by 
the  initial  condition  (8),  which  is: 


d£'fyx,y,0)  _  El  -Eq 
dt  At 


(16) 


Using  E]  and  Eq,  the  non- boundary  value  of  E2  is  given  by  equation  (14).  Then  the 
boundary  conditions  are  added  to  E2  and  the  process  is  repeated  to  find  E3  from  E2  and  Ei 
until  the  desired  time  is  reached. 
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Figure  7.  FDTD  space- time  node  diagram 
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B.  MODEL  DESCRIPTION 

During  the  forward  or  recording  step,  equation  (14)  is  used  to  find  the  forward 
time  solutions  of  the  wave  equation  using  the  initial  conditions  and  the  boundary 
conditions.  Equation  (14)  is  also  used  during  the  reverse  or  reconstruction  step  to  find 
the  reverse  time  solutions  of  the  wave  equation  by  defining  the  solutions  found  at  the 
desired  time  step  T  as  the  initial  and  the  boundary  conditions. 

Since  it  is  not  possible  to  consider  an  infinite  region  to  propagate  the  fields,  a 
limited  region  is  required  to  terminate  the  simulation.  Rather  than  utilizing  absorbing 
type  boundary  conditions  which  yield  only  approximate  cancellation  of  reflections,  it  was 
decided  to  avoid  all  reflections  by  employing  a  distant  terminating  boundary  whose 
reflection,  although  large,  does  not  arrive  back  at  the  sensor  boundary  until  t  >  T. 

1.  Model 

Denoting  1V|,  and  Np  as  the  number  of  receivers  on  the  recording  surface  of  the  x- 
and  y- segment,  and  Me  and  Ne  as  the  spacing  between  each  of  the  x-  and  y- segment 
receivers,  respectively,  the  number  of  grid  nodes  on  the  sensor  surface  segments  becomes 

Ms=(Mp- 1  )xMe+ 1  (x-  segment) 

Ns=(Np-l)xNe+l  (y- segment)  (17) 

The  number  of  nodes  on  the  distant  boundary  will  be  set  in  this  work  to  be 

M=2xMs  (x- segment) 

N=2xNs  (y- segment)  (18) 

The  sensor  surface  is  located  at  the  approximate  center  of  the  distant  boundary  by 
defining  its  edges  as  follows: 
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ml=  (M-Ms)/2+l  (Left  edge),  m2=  ml+Ms- 1  (Right  edge)  (x- segment) 


nl=  (N-Ns)/2+l  (Lower  edge),  n2=n2+Ns-l  (Upper  edge)  (y-segment)  (19) 
The  geometry  is  described  in  Figure  8. 

For  simphcity,  the  spacing  of  each  node  (h)  is  set  to  1  m  throughout  this  thesis.  As 
a  result,  the  2-D  stabihty  condition  (15)  is  rewritten  as 

^  >  ^/2  (20) 

cAt 


Usually,  in  two  or  three-  dimensional  problems  the  stabihty  condition  q  is  chosen 
to  be  twice  that  needed.  Therefore,  this  model  also  chooses  q  to  be  2.  This  ensures  two 
samples  between  two  adjacent  points  at  each  time  step. 

2.  Initial  And  Boundary  Conditions 

Since  this  thesis  considers  only  the  source- free  region,  no  fields  exist  within  the 
recording  surface  and  the  distant  boundary  at  the  initial  time.  Therefore,  the  initial 
conditions  within  the  distant  boundary  are  given  by: 


E^{m,n,l)  =  0 
dE^(m,n,l)  _ 


(21) 


where  m  =  \,M  and  n  =  1,  N. 

Assuming  the  scattered  field  is  to  be  reflected  at  the  distant  boundary,  the 
boundary  conditions  at  the  distant  boundary  become: 

E(l,n,p)  =  0,  E^(M,n,p)  =  0 

(22) 

E^(mXp)=0,  E^(m,N,p)  =  0 
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where  m  =1,  M,  n  =1,  N  and  p  =1,  T.  These  boundary  conditions  are  maintained 
throughout  the  forward  time  step. 

The  field  data  at  the  recording  surface  must  be  stored  at  each  time  step.  Letting 
EzBCl,  EzBC2,  EzBC3,  and  EzBC4  denote  the  arrays  storing  the  field  data  of  each  side  of 
the  recording  surface  during  the  forward  time  step: 

EzBC\{p)  =  E^{m^,n\,p),  EzBC2{p)=  E^(m^,n2,  p) 

EzBC3{p)  =  E  p),  EzBC4{p)=  E^(m2,n^,p) 

where  m^=  mi,  ml+Mg,  ml+2Me, . m2,  nr  =  nl,  nl+Ne,  nl+2Ne, . n2,  and p=l,T. 

The  field  data  stored  at  the  recording  surface  at  any  forward  time  step,  p  =  P  and  p 
=  P-1,  where  Pis  less  than  T,  can  be  defined  as  the  initial  conditions  for  the  time-reversed 
process: 


E[im^,nl,l)  =  EzBCliP),  E[{m^,n2,l)  =  EzBC2{P) 

E[{m\,n^  ,1)  =  EzBC3{P),  E[{m2,n^  ,1)  =  EzBC4(P) 

and 

E[{m^,nl,2)  =  EzBCl{P  -1),  E[{m^,n2,2)  =  EzBC2{P -1) 

E[  (ml,  n,  ,2)  =  EzBC3(P  - 1),  E[  (m2,  n,  ,2)  =  EzBC  4{P  - 1) 

where  E[  is  the  time-reversed  field  and  P  <  T.  Based  on  conditions  (24)  and  (25), 

equation  (14)  is  used  to  update  the  time -reversed  field  inside  the  recording  surface  until 
the  reversed  time  step,  P. 

If  IvT  =  Ne  =  1,  the  time-reversed  field  is  easily  updated  because  ?h  is  assumed  to 
be  1.  If  Me  ?  1  or  Ne  ?  I,  however,  the  field  data  requires  interpolation  at  grid  points 
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between  each  two  adjacent  transducers  to  maintain  the  continuity  of  the  field.  To  satisfy 
this  condition,  we  use  a  ‘cubic  spline’  interpolation  algorithm  in  MATLAB. 

3.  Incident  Field 

The  incident  field  used  here  is  represented  by  a  Gaussian  impulse  plane  wave 
propagating  in  the  -x  direction,  centered  at  x  =  at  t  =  0,  with  the  form, 

(26) 

where  c  is  the  standard  deviation  of  spatial  width  of  the  pulse.  Viewing  the  impulsive 
waveform  at  x  =  Xg  gives 

f(x„t)  =  e-^^‘'  (27) 


where 


2o 


.  Its  frequency  spectmm  is  given  by: 


f(^0’  /)  =  ^ 


(28) 


The  3dB  bandwidth  of  this  baseband  pulse  is  given  by  solution  of 


V2 


(29) 


which  yields. 


fidB  ~  _ 


A  ln2 


71  V  2  27ia 


(30) 


Note  the  reciprocal  relationship  between  the  spatial  standard  deviation  and  the  3dB 
bandwidth. 
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4.  Field  Representation 

A  scattered  field  results  when  an  incident  field  interacts  with  a  target.  Interactions 
include  polarization  currents  excited  within  dielectric  materials  and  surface  currents 
induced  near  to  the  surface  of  metahic  conductors.  Assuming  a  target  composed  of 
perfectly  conducting  metallic  surfaces,  the  total  tangential  E- field  is  forced  to  zero  at  the 
surface.  Since  the  total  field  is  given  ,  where  Etotai,  Emc,  and  Es 

represent  the  total  field,  the  incident  field  and  the  scattered  field,  respectively,  the 
tangential  scattered  field  is  equal  to  the  negative  of  the  incident  tangential  field  at  the 
metaUic  surface  of  the  target. 

To  properly  execute  the  FDTD  time -reversed  field  simulation  of  the  scattered 
field  when  the  field  focuses  on  a  metaUic  target  we  need  to  provide  the  negative  incident 
field  boundary  condition  at  all  surface  target  nodes.  However,  we  usually  do  not  know 
where  the  target  is  or  when  the  incident  field  reaches  the  target.  Unless  this  field 
boundary  condition  is  apphed  at  these  target  nodes  during  the  time-reversed  process,  field 
interactions  occur  between  adjacent  nodes  and  the  time-reversed  scattered  field  solution 
does  not  become  extinguished  for  times  prior  to  the  initial  plane  wave  impact. 

Suppose  that  a  point- like  target  exists  inside  the  recording  surface  and  the  total 
scattered  field  on  the  recording  surface  is  already  known.  The  total  scattered  field  inside 
the  recording  surface  would  then  consist  of  the  scattered  field  generated  by  the  time- 
reversed  boundary  conditions  on  the  recording  surface  and  from  the  negative  incident 
field  boundary  conditions  caused  by  nature  on  the  target,  or 

=  +  (31) 
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where  is  the  total  scattered  field  inside  the  surface,  is  the  scattered  field 

generated  by  the  recording  boundary  and  is  the  scattered  field  interaction  with  the 

target.  We  can  solve  for£|^  by  using  only  the  boundary  data  on  the  recording  surface 
without  knowing  the  location  of  the  target.  is  able  to  be  set  to  zero  before  the  time 

,to,  when  the  incident  field  reaches  the  target.  As  a  result,  we  can  predict  as: 

El-r{rj-t)=-El^{r^-t)  (32) 

where  iy  is  the  location  of  the  target.  Once  the  node  location  is  identified,  we  can 
simulate  the  target  node  interaction  and  properly  terminate  the  time-  reversed  process. 

5.  Conservation  Of  Energy 

The  above  sequential  node  identification  process  may  be  extended  in  theory  to 
multiple  target  nodes.  However,  a  different  concept  is  applied  in  this  thesis  to  yield  a  cost 
function  for  solution  optimization.  This  concept  is  that  of  conservation  of  energy.  The 
total  scattered  field  energy  accumulated  within  the  recording  surface  during  the  time- 
reversed  process  will  become  constant  at  all  times  prior  to  the  initial  impact  of  the 
incident  field  on  the  target.  This  same  result  will  occur  in  the  time-reversed  FDTD 
simulation  if  the  reversed  field  collapses  while  the  correct  target  node  locations  are  given 
negative  incident  field  boundary  conditions.  If  the  accumulated  total  energy  in  the  time- 
reversed  solution  continues  increasing  before  this  initial  impact  time,  then  one  or  more 
assumed  target  nodes  are  wrong. 

The  energy  metric  within  the  sensor  surface  at  an  arbitrary  time  step,  t,  is  defined 
by: 
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(33) 


m2  nl 

Energy{ti)=  (m,n,?i)]^ 

m=ml  n=nl 

where  E^(m  ,n  ,ti)  is  the  time-reversed  field  at  each  node  within  the  recording  surface  at 
the  time  step,  ti .  The  accumulated  total  energy  to  the  time  step,  T,  is  defined  by 

T  m2  n2 

Energy  {t)  -  (34) 

m=ml  n=nl 

If  the  accumulated  energy  in  (34)  becomes  constant,  this  imphes  that  the  time-reversed 
field  is  collapsing  onto  the  exact  target  node  locations. 

An  example  is  presented  in  the  next  chapter  to  illustrate  this  concept.  In  the 
example,  we  enforce  the  negative  incident  field  on  the  exact  target  nodes  during  the 
reversed  time  step  and  obtain  the  accumulated  total  energy.  Next,  the  assumed  target 
nodes  are  shghtly  perturbed  from  the  original  location  to  demonstrate  the  increase  of 
accumulated  total  energy. 
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Figure  8.  Simulation  model  diagram 
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IV. 


EXAMPLES  OE  TIME-REVERSED  IMAGING 


This  chapter  presents  fundamental  examples  of  the  time -reversed  imaging 
process.  The  first  example  illustrates  the  basic  time-reversed  process,  while  the  other 
examples  consider  the  process  working  under  practical  consideration,  such  as  pulsed 
radar  applications.  However,  this  chapter  does  not  consider  some  aspects  of  the  real 
environment,  such  as  noise,  actual  3-D  shapes  and  materials  of  targets  and  any  kind  of 
information  loss  by  the  media.  In  aU  examples,  the  dimensions  of  the  FDTD  grid  and 
number  of  time  steps  are  set  as  follows:  IVfe  =  Ns=  101,  T  =  300.  This  selection  sets  the 
recording  surface  as  a  100m  x  100m  square  grid  and  the  distant  boundary  as  a  200m  x 
200m  square  grid  with  201  by  201  points. 

AU  programs  used  for  calculations  shown  in  this  chapter  are  presented  in  the 
Appendix. 

A.  TARGETS  NODES  WITH  INITIAL  CONDITIONS  ONLY 

This  first  example  considers  two  aircraft-like  targets  present  inside  the  recording 

surface.  The  targets  are  represented  by  specific  grid  points  that  are  excited  by  initial 
conditions,  Ez  (x,  y,  t=0)  =1,  but  with  no  subsequent  boundary  conditions.  Grid  points 
surrounding  the  target  nodes  are  also  endowed  with  initial  conditions,  but  which  rapidly 
decrease  with  distance,  per  a  Gaussian  distribution.  No  incident  field  is  considered  and 
there  are  no  boundary  field  interactions  for  the  target  nodes,  such  as  with  conducting 
targets.  This  example  demonstrates  how  the  reversed  field  converges  onto  the  target 
nodes  for  the  fundamental  initial  value  problem. 
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1.  Forward  Time  Step  Solutions 

The  forward  time  simulation  is  described  first  to  help  understand  the  reversed 
time  simulation  described  later.  Figure  9  shows  the  initial  condition  of  the  simulation. 
For  convenience,  the  recording  surface  is  presented  as  a  dashed  hne.  Every  60th  time  step 
is  presented  in  Figure  10.  Figure  11  represents  the  final  time  step.  Time-reversed 
recorded  field  values  on  the  inner  grid  boundary  are  used  to  drive  the  time-reversed 
solution  while  null  initial  conditions  are  assumed  within  the  inner  grid. 


Log.,Q{|E^(x,y.l=1)|};  Fofward  Time  Steps:  300 


x-axis 


Figure  9.  Initial  condition  of  the  forward  time  step  for  two  aircraft-tike  targets 
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Figure  10.  Forward  time  step  for  two  aircraft -like  targets  atT  =  61,T  =  121, T  =  181  and 
T  =  241 
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Figure  1 1 .  Final  time  step  T  =  300  for  two  aircraft- like  targets 
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2.  Reversed  Time  Step  Solution  with  Me  =  Ne  =  1 

This  subsection  presents  the  reversed  time  step  simulation  in  the  case  of  IVfc  =  Ne  = 
1.  In  this  case,  each  node  acts  as  a  transducer;  thus,  the  number  of  transducers  is  Mp  =  Np 
=  101.  Selecting  Me  =  Ne  =  1  provides  the  most  accurate  field  representation  although 
some  errors  occur  due  to  the  computational  truncation.  Therefore,  the  final  reversed  time 
step  condition  does  not  indicate  the  exact  initial  condition  of  the  forward  time  step.  As 
seen  in  Figure  11,  the  field  reflected  from  the  outer  boundary  has  crossed  the  inner 
recording  surface  sometime  before  the  T=300  time  step.  This  reflection  must  not  be 
used  in  the  reverse-time  process.  Thus,  the  reverse-time  processing  will  begin  at  T=250. 
Figure  12  shows  the  initial  condition  of  the  reversed  time  step  with  T  =  250.  Every  25th 
reversed  time  step  is  displayed  in  Figure  13  and  Figure  14,  with  T  =  25  and  T  =  0  shown 
in  Figure  15. 

Log ^g{|E^(x,y, 1*250)1}  Initial  Condition  for  Time-Reversal 


Figure  12.  Initial  condition  of  the  reversed  time  step  for  two  aircraft- tike  targets  with 
Me=  Ne=  1  case 
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Figure  13.  Reversed  time  step  for  two  aircraft-like  targets  with  Me=  Ne=  1  at  T=  225,  T= 
200,  T=  175  and  T=  150 
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Figure  14.  Reversed  time  step  for  two  aircraft- like  targets  with  Me  =  Ne  =  1  at  T=  125, 
T=  100,  T=  75  and  T=  50 
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Figure  15.  Reversed  time  step  for  two  aircraft- tike  targets  with  Me  =  Ne=  1  at  T=  25  and 
T=0 
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3.  Reversed  Time  Step  Solution  with  Me  =  Ne  =  5 

This  subsection  presents  the  case  of  Me  =  Ne  =  5,  which  corresponds  to  Mp  =  Np= 
21.  This  selection  means  that  every  5th  point  data  on  the  recording  surface  is  used  to 
calculate  the  reversed  field.  Figure  16  shows  the  initial  condition  of  the  reversed  time 
step  at  T  =  250.  The  red  circles  in  Figure  16  represent  the  transducers.  Figure  17  and  18 
display  every  25th  time  step,  while  Figure  19  shows  the  final  condition  of  the  reversed 
time  step  at  T=0.  This  Me  =  Ne=  5  case  provides  less  accuracy  than  the  case  of  Me  =  Ne  =1 
does,  but  it  still  gives  useful  targets  information. 


Log ^Q{|E^(x.y, 1*250)1}  Initial  Condition  for  Time-Reversal 
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Figure  16.  Initial  condition  of  the  reversed  time  step  for  two  aircraft- tike  targets  with  Me 
=  Ne=5 
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Figure  17.  Reversed  time  step  for  two  aircraft -like  targets  with  Mg  =  Ne  =  5  at  T=  225, 
T=  200,  T=  175  and  T=  150 
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Figure  18.  Reversed  time  step  for  two  aircraft- like  targets  with  Me  =  Ne  =  5  at  T=  125,  T= 
100,  T=  75  and  T=  50 
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Figure  19.  Reversed  time  step  for  two  aircraft- like  targets  with  Me  =  Ne  =  5  at  T=  25  and 
T=0 
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4.  Reversed  Time  Step  Solution  With  Me  =  Ne=  10 

In  this  final  section,  the  case  of  Me  =  Ne=  10  case  is  simulated.  In  this  situation, 


the  corresponding  number  of  the  transducers  is  Mp  =Np  =11  indicating  that  every  10th 
point  on  the  recording  surface  is  used.  Figure  20  shows  the  initial  condition  of  the 
reversed  time  step  at  T  =  250.  Figure  21  and  22  display  every  25th  time  step.  Figure  23 
shows  the  final  step  of  the  reversed  time  step,  which  makes  the  targets’  shape  difficult  to 
identify.  In  this  case,  the  images  of  the  targets  are  blurred. 

±  Log^Q{|E^(x.y.t*250)|}  Initial  Condition  for  Time-Reversal 


Figure  20.  Initial  condition  of  the  reversed  time  step  for  two  aircraft- like  targets  with 

Me=  Ne=  10 
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Figure  21 .  Reversed  time  step  for  two  aircraft- like  targets  with  Me  =  Ne  =  10  at  T=  125, 
T=  100,  T=  75  and  T=  50 
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Figure  22.  Reversed  time  step  for  two  aircraft- like  targets  with  Me  =  Ne=  10  at  T=  125, 
T=  100, T=  75  and  T=  50 
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Figure  23.  Reversed  time  step  for  two  aircraft- like  targets  with  Me  =  Ne  =  10  at  T=  25 
and  T=  0 
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B.  TWO  POINT-LIKE  SCATTERING  TARGETS 

In  the  previous  example  of  two  aircraft- like  targets  an  outbound  field  is  generated 

due  only  to  initial  conditions  without  an  incident  field.  In  this  section  we  consider  an 
impulse  radar  example  with  scattering  of  a  Gaussian  plane  wave  by  two  point- like 
metaUic  targets.  The  Gaussian  plane  wave  employed  in  this  model  has  a  standard 
deviation,  s  =lm,  which  corresponds  to  a  3dB  effective  bandwidth,  BWsdB  ~  125  MHz. 
In  this  example,  two  point- like  targets  are  used  to  process  the  reversed  time  operation. 
Again,  the  size  of  the  recording  surface  is  IW  =  Ns  =  101,  which  is  a  100m  x  100m  square 
grid  for  the  recording  surface  and  a  200m  x  200m  square  grid  for  he  distant  boundary. 
The  forward  time  step  solutions  are  presented  in  the  next  subsection  followed  by  the 
reversed  time  step  solutions. 

I.  Forward  Time  Step  Solutions  For  Two  Point -Like  Targets 

This  subsection  presents  the  forward  time  step  solution  generated  by  an  impulsive 

plane  wave  incident  field.  Boundary  data  at  the  recording  surface  is  recorded  at  each 
time  step.  Two  point- hke  perfect  conductor  targets  forcing  the  total  B  field  to  $ro  exist 
at  the  approximate  center  of  the  grid.  The  incident  field  starts  marching  from  x  =120. 
For  convenience,  the  targets  are  represented  by  two  black  dots  in  all  figures  and  the 
recording  surface  is  displayed  as  dashed  lines.  Figure  24  shows  the  initial  condition  of 
the  forward  time  step  solution  for  two  point- hke  targets  with  an  incident  field.  Figure25 
displays  the  results  at  T=  51,T=  101,T=  151  and  T=  201.  Finally,  the  result  of  T=  251  and 
the  final  condition  at  T=  300  are  presented  on  Figure  26. 
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Figure  24.  Initial  condition  of  the  forward  time  step  solution  for  two  point-like  targets 
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Figure  25.  Forward  time  step  for  two  point- like  targets  at  T=  51,  T=  101,  T=  151  and  T= 
201 
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Figure  26.  Forward  time  step  for  two  point- like  targets  at  T=  25 1  and  T=  300 
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2.  Reversed  Time  Step  Solutions  For  The  Two  Point -Like  Targets 

This  subsection  provides  the  results  of  the  reversed  time  step  simulation.  To 

simulate  the  exact  time-reversed  solution  we  must  enforce  the  PEC  boundary  condition  of 
zero  total  Bfield  at  each  target  node.  Our  goal  is  to  recover  the  targets’  shape,  which  is 
assumed  to  be  unknown,  so  determining  the  correct  nodes  to  enforce  the  PEC  boundary 
condition  will  be  akin  to  solving  the  imaging  problem.  A  correct  selection  of  the  target 
nodes  will  make  the  reversed  scattered  field  collapse  onto  the  targets  and  become 
completely  extinguished  at  all  times  prior  to  the  initial  impact  of  the  incident  plane  wave. 

As  seen  in  the  forward  time  step  results,  the  time  step,  T  =  200,  is  a  good  choice 
to  begin  the  reversed- time  simulation.  In  this  example  the  spacing  of  each  transducer  is 
assumed  to  be  1:  Me  =  Ne  =  1.  Only  the  scattered  field  is  displayed  in  this  example.  Eor 
convenience  a  red  dash  tine  is  used  to  represent  the  location  of  the  incident  field  peak  at 
each  time  step.  The  two  point- like  targets,  represented  by  two  black  dots,  are  located  at 
Target  1:  (x,  y)  =  (96,101)  and  Target  2:  (x,  y)  =  (106,101).  During  the  simulation,  the 
scattered  field  is  forced  to  be  the  negative  incident  field  at  the  target  nodes.  Because  the 
simulation  uses  the  data  as  exactly  recorded  on  the  recording  surface  during  the 
corresponding  forward  time  step,  the  scattered  field  collapses  on  the  targets.  Eigure  27 
shows  the  initial  condition  of  the  reversed  time  step  at  T  =  200.  Eigure  28  and  29  display 
every  25th  time  step.  Eigure  30  indicates  the  accumulated  total  energy  inside  the 
recording  surface,  where  the  xaxis  corresponds  to  the  reversed  time  step,  T  =  200  -  twith 
t  being  the  forward  time  step.  The  accumulated  total  energy  inside  the  recording  surface 
remains  constant  after  the  scattered  field  collapses  onto  the  targets. 
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Figure  27.  Initial  condition  of  the  reversed  time  step  for  the  exact  two  point- tike  targets 
at  T  =  200 
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Figure  28.  Reversed  time  step  for  the  exact  two  point-like  targets  at  T  =  175,  T  =  150, T 
=  125  andT=  100 
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Figure  29.  Reversed  time  step  for  the  exact  two  point-like  targets  at  T  =  75,  T  =  50, T  = 
25  and  T  =  1 
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Final  Sub-Grid  Energy  =  1346.7036 


Figure  30.  Accumulated  total  energy  inside  the  recording  surface 
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3.  Reversed  Time  Step  Solutions  For  Perturbed  Two  Point -like  Targets 

This  subsection  provides  the  reversed  time  solutions  for  PEC  boundary  conditions 

enforced  at  locations  slightly  in  error  from  the  exact  nodal  locations  of  the  two  point- like 
targets. 

The  negative  incident  field  is  enforced  at  nodes  represented  by  red  dots  while  the 
exact  target  node  locations  are  again  indicated  by  black  dots.  Since  the  data  used  for  the 
reversed  time  step  is  based  on  the  forward  time  step  solutions  in  which  the  targets  are 
located  at  the  correct  nodes,  the  time-reversed  scattered  field  does  not  become 
extinguished,  but  diverges  after  the  reversed  field  reaches  the  new  targets.  Furthermore, 
the  accumulated  total  energy  inside  the  recording  surface  keeps  increasing  although  the 
incident  field  initially  does  not  reach  the  targets  in  the  forward  time  step.  To  confirm  the 
concept  of  the  accumulated  total  energy,  three  cases  are  examined.  In  the  first  case. 
Target  1  is  perturbed,  but  Target  2  is  not  perturbed.  In  the  second  case.  Target  I  stays  at 
the  original  coordinate,  while  Target  2  is  perturbed.  In  the  last  case,  both  targets  are 
perturbed. 

a.  Case  1:  Target  1  Is  Perturbed,  Target  2  Is  Not 

In  this  case.  Target  I  is  perturbed  to  (x,  y)  =  (90,101),  but  Target  2  is 

stationed  at  the  correct  location,  (x,  y)  =  (106,  101).  The  initial  condition  at  T  =  200  is 

shown  in  Figure  31.  Figure  32  displays  the  results  of  every  40th  time  step.  Figure  33 

shows  the  final  condition  at  T  =  1.  The  accumulated  total  energy  inside  the  recording 

surface  is  shown  in  Figure  34.  This  case  causes  the  accumulated  total  energy  to  increase 

after  the  reversed  field  reaches  the  new  targets. 
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b.  Case  2:  Target  1  Is  Not  Perturbed,  Target  2  Is  Perturbed 

In  this  case,  Target  1  is  located  at  the  correct  coordinate,  (x,  y)  =  (96,  101), 

but  Target  2  is  perturbed  to  (x,  y)  =  (110,101).  The  initial  condition  at  T=200  is  shown  in 
Figure  35.  Figure  36  displays  every  40th  time  step  and  Figure  37  shows  the  final 
condition  at  T=l.  Figure  38  shows  the  accumulated  total  energy. 

c.  Case  3:  Both  Targets  Are  Perturbed 

Target  1  is  perturbed  to  (x,  y)  =  (90,101)  and  Target  2  is  perturbed  to  (x, 
y)=(110,  101).  Figure  39  shows  the  initial  condition  at  T  =  200.  Figure  40  displays  every 
40th  time  step  and  Figure  41  shows  the  final  condition  at  T  =  1.  The  accumulated  total 
energy  is  plotted  in  Figure  42. 
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Figure  3 1 .  Initial  condition  of  the  reversed  time  step  for  Case  1 
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Figure  32.  Reversed  time  step  solutions  for  Case  1  at  T  =  160,  T  =  120,  T  =  80  and  T  = 
40 
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Figure  33.  The  final  condition  of  the  reversed  time  step  for  Case  1 
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Final  Sub-Grid  Energy  =  1619.2719 


Figure  34.  Accumulated  total  energy  inside  the  recording  surface  for  Case  1 
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Figure  35.  Initial  condition  of  the  reversed  time  step  for  Case  2 


58 


30Q 


•  LogtollE^^.y.t^tO)!) 


x-axIs 


Figure  36.  Reversed  time  step  solutions  for  Case  2  at  T  =  160,  T  =  120,  T  =  80  and  T  = 
40 
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Figure  37.  The  final  condition  of  the  reversed  time  step  for  Case  2 
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Final  Sub-Grid  Energy  =  1483.6283 


Figure  38.  Accumulated  total  energy  inside  the  recording  surface  for  Case  2 


61 


±  Log,Q{|E^«'(x.y.l“200)|} 


Figure  39.  Initial  condition  of  the  reversed  time  step  for  Case  3 
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Figure  40.  Reversed  time  step  solutions  for  Case  3  at  T  =  160,  T  =  120,  T  =  80  and  T  = 
40 
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Figure  41 .  The  final  condition  of  the  reversed  time  step  for  Case  3 
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Final  Sub-Grid  Energy  =  1774.035 


Figure  42.  Accumulated  total  energy  inside  the  recording  surface  for  Case  3 
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C.  SUMMARY 

This  chapter  presented  the  fundamental  examples  of  the  time-reversed  process  for 
electromagnetic  waves  by  using  the  FDTD  method.  The  first  example  demonstrated  the 
time-reversed  process  when  the  target  nodes  radiate  due  only  to  enforced  initial 
conditions  without  an  incident  field.  Essentially  perfect  reconstmction  of  the  two 
aircraft- tike  target  images  is  obtained  when  all  nodes  on  the  recording  surface  are  used  to 
drive  the  time-reversed  solution.  The  effect  of  thinned  spacing  between  the  transducers 
was  then  examined.  This  example,  however,  is  not  related  to  the  radar  scattering  case. 

The  next  example,  which  considered  an  incident  Gaussian  impulse  plane  wave, 
showed  that  the  time-reversed  scattered  field  collapses  onto  the  point-tike  targets,  with 
complete  extinction  only  if  the  PEC  boundary  condition  is  enforced  on  the  target  nodes 
during  the  reversed  time  step.  A  wrong  selection  for  any  target  node  location  yields  a 
reversed- time  scattered  field  that  initially  collapses  but  then  expands  outward.  Incorrect 
target  node  selections  provide  scattered  field  energy  at  times  earlier  than  the  initial 
impact  of  the  incident  field.  The  accumulated  energy  in  the  scattered  field  continues  to 
increase  for  any  wrong  target  node  selection  as  shown  in  the  three  cases  examined.  It 
thus  appears  that  accumulated  energy  in  the  reversed- time  scattering  solution  can  be 
employed  as  a  “cost  function”  in  an  optimization  routine  designed  to  find  the  correct 
target  nodes. 
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V.  SUMMARY  AND  CONCLUSIONS 


A.  SUMMARY 

This  thesis  presents  an  initial  attempt  to  demonstrate  target  imaging  by  using  a 
numerical  time-reversed  process  for  the  radar  imaging  applications,  in  particular,  for  bi¬ 
static  or  multi- static  radars. 

Demonstration  of  the  time -reversed  process  is  done  by  the  FDTD  method.  The 
first  example  shows  that  the  time-reversed  process  can  image  multiple  targets  when  the 
only  driving  function  for  the  field  is  initial  conditions  at  the  target  nodes.  This,  however, 
is  not  the  radar  scattering  case. 

The  radar  scattering  case  was  examined  next,  albeit  for  a  simple  example  of  two 
point  targets  illuminated  by  an  incident  Gaussian  pulse  plane  wave.  A  perfect  time- 
reversed  simulation  provides  a  scattered  field  that  collapses  onto  the  target  nodes  and 
becomes  extinguished  for  aU  times  earlier  than  the  initial  plane  wave  impact.  Such  a 
solution  requires  enforcement  of  metallic  boundary  conditions  (scattered  field  equals 
negative  of  the  incident  field)  at  each  target  node.  This,  of  course,  requires  a  priori 
knowledge  of  the  target  location,  orientation  and  geometry  -  the  very  knowledge  being 
sought  in  the  time-reversed  process.  AU  is  not  lost  however,  by  what  appears  to  be  a 
chicken  and  egg  dilemma.  Wrong  selections  of  even  one  target  node  for  enforcing  the 
metaUic  boundary  condition  yields  time-reversed  scattering  solutions  which  first  coUapse 
then  radiate  outward  from  the  target  region  rather  than  being  extinguished.  The 
accumulated  energy  in  the  solution  can  thus  be  used  as  a  cost  function  for  optimization  of 
the  scattering  nodes,  with  the  minimum  cost  solution  corresponding  to  the  correct 
locations. 
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B.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  function  of  a  radar  imaging  system  is  to  gain  early  identification  of  non- 

cooperative  targets.  The  objective  of  this  thesis  was  to  initially  investigate  the  potential 
utility  of  the  time- reversed  numerical  process  for  employment  in  multi- static  radar 
imaging  applications.  Using  an  initial  condition  solution  as  the  first  case  to  investigate, 
it  was  found  that  superb  imaging  is  possible  for  multiple  targets  when  all  boundary 
nodes,  without  thinning,  are  used  to  drive  the  time-reversed  solution.  The  radar 
scattering  case,  however,  is  a  different  animal  that  will  require  solution  of  an 

optimization  problem  where  the  target  nodes  are  the  parameters  to  be  optimized 
(correctly  placed)  and  the  cost  function  is  the  accumulated  scattered  field  energy  in  the 
time-reversed  numerical  solution. 

Future  work  on  this  subject  should  include  investigation  of  optimization 

algorithms  for  solution  of  the  correct  target  nodes.  The  genetic  algorithm  may  be  the 

best  bet  for  this  task  since  it  is  robust  in  exploring  a  wide  range  of  local  cost  minima  in 
solving  for  the  best  global  solution  within  the  parameter  space.  In  addition,  a  three- 
spatial  dimension  model  supporting  volumetric  regions  should  be  developed.  To 

represent  a  reahstic  environment,  this  model  should  include  provision  for  uneven  spacing 
of  the  transducers  with  3-D  targets.  The  model  should  be  evaluated  relative  to  required 
operating  frequencies,  signal- to- noise  ratio  and  the  number  and  material  properties  of 
realistic  targets. 
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APPENDIX.  PROGRAM  LISTINGS 


A.  PROGRAM  LISTINGS 

These  programs  are  used  to  demonstrate  the  forward  scattering  and  the  time- 
reversed  process  in  MATLAB.  The  program  used  for  the  first  example  allows  inputs  for 
the  number  of  transducers,  the  spacing  between  the  transducers,  the  stabihty  condition, 
time  step,  the  initial  condition  of  the  region  and  the  dynamic  range  for  the  plot.  The 
program  used  for  the  second  example  allows  inputs  for  the  number  of  transducers,  time 
step,  the  initial  position  of  the  incident  field  and  the  standard  deviation  of  the  incident 
field.  Target  functions  are  also  hsted. 


These  programs  are  written  by  Prof.  M.  A.  Morgan  and  modified  by  LT.  Inaba. 

I.  Program  RTWE_2D5  For  The  First  Example 

Mp=mput('Enter  number  of  field  sensors  on  each  x-segment  of  the  sub-grid  boundary:  '); 
Me=input('Enter  node  spacing  between  each  x-segment  subgrid  boundary  sensor: '); 
Np=input('Enter  number  of  field  sensors  on  each  y-segment  of  the  sub-grid  boundary: '); 
Ne=input('Enter  node  spacing  between  each  y-segment  subgrid  boundary  sensor:  '); 

Ms=(Mp-l)*Me  -r  1;  Ns=(Np-l)*Ne  -r  1; 
dispd'Subgrid:  Ms=',mt2str(Ms),';  Ns=',int2str(Ns)]); 

%  Default  Mesh 
M=2*Ms;  N=2*Ns; 


ml=fix((M-Ms)/2)-i-l;  m2=ml-i-Ms-l; 
nl=fix((N-Ns)/2)-rl;  n2=nl-rNs-l; 
x=l:M;  y=l:N;  Y=y'*ones(l,M);  X=ones(N,l)*x; 
xs=ml:m2;  ys=nl:n2;  Ys=ys'*ones(l,Ms);  Xs=ones(Ns,l)*xs; 
dispd'Grid:  M=',mt2str(M),'  ml=',int2str(ml),'  m2=',int2str(m2),... 

N=',int2str(N),'  nl=',mt2str(nl),'  n2=',mt2str(n2)]); 
xe=(ml:Me:m2)';  ye=(nl:Ne:n2)';  %  sub-grid  sensor  x  and  y  node  numbers 

dispCNote  q  >=  sqrt(2)  Courant  requirement  for  convergence  '); 
q=input('Enter  Integer  Value  for  q=dh/(c*dt)  (1,  2,  etc):  '); 

Nmin=min(  1 .5  *M-m2, 1 .5  *N-n2) ;  Pmax=fix(q*Nmin) ; 

P=inputd'Enter  number  of  time-steps  (  <  ',int2str(Pmax),')  :  ']); 
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namel=input('Enter  Forward  Soln  Files  "Name"  for  Name_nn.jpg  ?  (Enter  Key  to  Skip): 
if  ~isempty(namel), 

dpl=input('Enter  Time  Step  Increment  Between  Stored  Frames:  '); 
pjl=l:dpl:P; 
nj= 1 ;  njmx=length(pj  1 ) ; 
end 

Ql=l/(q*q);  Q2=(2-4*Q1); 

Ez=zeros(N,M,P);  %  Reserving  3D  array 
EzIC=zeros(N,M) ; 

EzBCl=zeros(Ms,P);  EzBC2=zeros(Ns,P);  EzBC3=EzBCl;  EzBC4=EzBC2; 

%  Subgrid  BC's 

%  Constmcting  Forward-Time  Evolution  Arrays 
d=ones(  1  ,max(M,N)) ; 

A=diag(d(l  :N-  l),l)+diag(d(l  :N- 1),- 1); 

B=diag(d(  1  :M  - 1 ),  1  )+diag(d(  1  :M  - 1 ) ,- 1 ) ; 

%  Define  Basic  Aircraft  Shaped  Grid  Focations  Centered  at  y=x=0 
[ny  ,mx]=Air2(N,M) ; 

Ntgt=length(ny); 

%  Defining  IC  Node  Centers 
for  nic=l:Ntgt 
EzIC(ny(nic),mx(nic))=l ; 
end 

%  Define  Unit  Peak  Gaussian  Shaped  Initial  Condition 
s=input('Enter  IC  Std  Dev  in  Grid  Spaces  (0  to  use  point  ICs):  '); 
if  ~isempty(s), 
if  s  >  0, 

[ny  mx]=fmd(EzIC==l);  %  Array  locations  for  IC  nodes 

Nic=length(ny); 

for  n=l:Nic 

Xc=X-x(mx(n));  Yc=Y-y(ny(n)); 

R2=(l/(2*s*s))*(Xc.*Xc-tYc.*Yc); 

Ez(:,:,l)=Ez(:,:,l)  -i-  exp(-R2);  %  Gaussian  spread  about  each  IC  node 
end 
else, 

Ez(:,:,l)=EzIC;  %  Using  unit  node  ICs  for  s  <=  0 

end 
end 

if  isempty(s),  Ez(:,:,l)=EzIC;  end 

%  Initial  Subgrid  BC's 
ifMe==  1, 

EzBCl(:,l)=Ez(nl,xs,l).’;EzBC3(:,l)=Ez(n2,xs,l).’; 

else, 

EzBCl  (:,  1  )=spline(xe,Ez(nl  ,xe,  l).',xs'); 

EzBC3(:,l)=spline(xe,Ez(n2,xe,l).',xs'); 
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end 

ifNe==  1, 

EzBC2(: ,  1  )=Ez(ys,m2, 1) ;  EzBC4(:, l)=Ez(ys,ml ,  1 ); 
else, 

EzBC2(: ,  1  )=spline(ye,Ez(ye,m2, 1  ),ys'); 

EzBC4(:,  1  )=spline(ye,Ez(ye,ml ,  l),ys'); 
end 

%  Assuming  dEz/dt=G=0  at  p=l  to  take  initial  step  to  p=2  (see  5/10/00  note) 
Ez(:,:,2)=0.5*(Ql*(A*Ez(:,:,l)  +  Ez(:,:,l)*B)  +  Q2*Ez(:,:,l)); 

%  p=2  Subgrid  BC's 
ifMe==  1, 

EzBCl(:,2)=Ez(nl,xs,2).’;EzBC3(:,2)=Ez(n2,xs,2).’; 

else, 

EzBC  1  (:  ,2)=spline(xe,Ez(nl  ,xe,2).',xs') ; 
EzBC3(:,2)=spline(xe,Ez(n2,xe,2).',xs'); 
end 

ifNe==  1, 

EzBC2(:,2)=Ez(ys,m2,2);  EzBC4(:,2)=Ez(ys,ml,2); 
else, 

EzBC2(:,2)=spline(ye,Ez(ye,m2,2),ys'); 

EzBC4(:,2)=spline(ye,Ez(ye,ml,2),ys'); 

end 

for  p=3:P;  %  Equation  of  Evolution 

Ez(:,:,p)=Ql*(A*Ez(:,:,p-l)+Ez(:,:,p-l)*B)+Q2*(Ez(:,:,p-l)-Ez(:,:,p-2)); 

%  Explicitely  Enforce  Ez=0  BC's  on  Grid  Boundary  Eor  Update 
Ez(  1 , :  ,p)=zeros(  1  ,M) ;  Ez(N, :  ,p)=zeros(  1  ,M) ; 

Ez(:,l,p)=zeros(N,l);  Ez(:,M,p)=zeros(N,l); 

Ezmx(p)=max(max(abs(Ez(: , :  ,p)))) ; 

%  Subgrid  BC  Update 
ifMe==  1, 

EzBC  1  ( :  ,p)=Ez(n  1  ,xs ,p). ' ;  EzBC3 (:  ,p)=Ez(n2,xs,p) ; 
else, 

EzBC  1  (:  ,p)=spline(xe,Ez(nl  ,xe  ,p) . '  ,xs') ; 
EzBC3(:,p)=spline(xe,Ez(n2,xe,p).',xs'); 
end 

ifNe==  1, 

EzBC2(:  ,p)=Ez(ys,m2,p) ;  EzBC4( :  ,p)=Ez(ys,m  1  ,p) ; 
else, 

EzBC2(:,p)=spline(ye,Ez(ye,m2,p),ys'); 

EzBC4( :  ,p)=spline(ye,Ez(ye,m  1  ,p),ys') ; 
end 
end 

while  1 

dispCMovie  Modes: ') 
disp('l  ->  2-D  log-scale  color  plan  view') 
disp('2  ->  2-D  linear  scale  color  plan  view') 
disp('3  ->  3-D  linear  scale  copper  tint') 
vsel=input('Select  Choice:  '); 
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if  ~isempty(vsel),  if  vsel==l  I  vsel==2  I  vsel==3,  break;  end;  end 
end 

EzMax=niax(Ezmx) ; 


if  vsel==l, 

%  Using  bi-polar  log  scaling  to  retain  Ez  polarity  with  selected  dynamic  range 
DR=input('Enter  2-D  Plot  Dynamic  Range,  e.g.  100,  1000,...  :  '); 
SEac=EzMax/DR;  Cmax=logl0(DR);  C=[-Cmax  Cmax];  Zc=2*Cmax; 
XBC=[ml  m2  m2  ml  ml]’;  YBC=[nl  nl  n2  n2  nl]’;  ZBC=Zc*[l  1  1  1  1]’; 

%  sub-grid  outline 

XBS=[xe;  xe;  ml*ones(Np,l);  m2*ones(Np,l)]; 

%  sub-grid  sensor  x  and  y  node  numbers 
YBS=[nl*ones(Mp,l);  n2*ones(Mp,l);  ye;  ye]; 

ZBS=Zc*ones(2*(Mp-tNp),  1 ); 
end 


if  vsel==2, 

%  Using  Linear  Ez  plot 
C=[-EzMax  EzMax];  Zc=2*EzMax; 

XBC=[ml  m2  m2  ml  ml]’;  YBC=[nl  nl  n2  n2  nl]’;  ZBC=Zc*[l  1  1  1  1]’; 
%  sub-grid  outline 

XBS=[xe;  xe;  ml*ones(Np,l);  m2*ones(Np,l)]; 

%  sub-grid  sensor  x  and  y  node  numbers 
YBS=[nl*ones(Mp,l);  n2*ones(Mp,l);  ye;  ye]; 
ZBS=Zc*ones(2*(Mp-tNp),  1) ; 
end 

if  vsel==3,  v(5)=-EzMax;  v(6)=EzMax;  end 
%  Initializing  Erames 

elf  reset;  v(l)=l;  v(2)=M;  v(3)=l;  v(4)=N; 


ifvsel==l, 

EzScl=Ez(:,:,l)/SEac;  %  scaling  so  abs(EzScl)  <=  DR 

[pi  ql]=fmd(0  <=  EzScl  &  EzScl  <  1);  %  nonlinear  remapping  for  lEzScll  <  1 
EzScl(pl-i-(ql-l)*N)=  1;  %  to  zero  log  plot 

[pi  ql]=fmd(-l  <  EzScl  &  EzScl  <  0); 

EzScl(pl+(ql-l)*N)=-l; 

EzLog=sign(Ez(: , : ,  1 )).  *log  1 0(abs(EzScl)) ; 

surf(X,Y,EzLog);  shading  interp 

title( [’Initial  Condition  \pm  Log_{  10}\[  IE_z(x,y,t=0)l\};  Eorward  Time  Steps:  ’... 
,int2str(P)]  ,’EontSize’ ,  1 4) 

xlabel(’x-axis’,’PontSize’,  1 8);  ylabel(’y-axis’,’PontSize’,  1 8) 

axis  equal;  axis(v);  caxis(C);  colorbar 

view(0,90);  hold  on 

plot3(XBC,YBC,ZBC,’  -k’) ; 

if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,’or’);  end 

hold  off 
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end 


if  vsel==2, 

surf(X,Y,Ez( shading  interp 

tide( ['Initial  Condition  E_z(x,y,t=0);  Eorward  Time  Steps: 

,int2str(P)]  ,'EontSize' ,  1 4) 

xlabel('x-axis','EontSize',  1 8);  ylabel('y-axis','EontSize',  1 8) 
axis  equal;  axis(v);  caxis(C);  colorbar 
view(0,90);  hold  on 
plot3(XBC,YBC,ZBC,’ -k’) ; 
if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,’or’);  end 
hold  off 
end 

if  vsel==3,  surfl(Y,X,Ez(:,:,l)); 

title( ['Initial  Condition  E_z(x,y,t=0);  Eorward  Time  Steps:  ',int2str(P)],... 
'PontSize',14) 

xlabel('x-axis','PontSize',  1 8);  ylabel('y-axis','PontSize',  1 8) 
axis(v);  colormap(copper); 
end 

Prame=moviein(P,gcf) ; 

hcpy=input('Print  a  Hard  Copy  ?  (Y/N):  ','s'); 
if  ~isempty(hcpy),  if  hcpy  ==  'Y'  I  hcpy  ==  'y',  print;  end;  end 
name=input('Enter  "Name"  to  Save  as  Name.jpg  Pile  ?  (Hit  Enter  Key  to  Skip):  ','s') 
if  ~isempty(name),  eval(['print  ',name,'  -djpeg99']);  end 

for  p=l:P;  %  Recording  animation  frames 
if  vsel==l, 

EzScl=Ez(:,:,p)/SPac;  %  scaling  so  abs(EzScl)  <=  DR 

[pi  ql]=fmd(0  <=  EzScl  &  EzScl  <  1);  %  nonlinear  remapping  for  lEzScll  <  1 
EzScl(pl+(ql-l)*N)=  1;  %  to  zero  log  plot 

[pi  ql]=find(-l  <  EzScl  &  EzScl  <  0); 

EzScl(pl+(ql-l)*N)=-l; 

EzLog=sign(Ez(: , :  ,p)).  *log  1 0(abs(EzScl)) ; 
surf(X,Y,EzLog);  shading  interp 

title(['\pm  Log_{  10}\{IE_z(x,y,t=',int2sfr(p),')l\};  Eorward  Time  Steps: 
int2str(P)],'PontSize',  14) 

xlabel('x-axis','PontSize',  1 8);  ylabel('y-axis','PontSize',  1 8) 

axis  equal;  axis(v);  caxis(C);  colorbar 

view(0,90);  hold  on 

plot3(XBC,YBC,ZBC,' -k') ; 

if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,'or');  end 

hold  off 

if  ~isempty(namel), 

if  p  ==  pjl(nj),  name=[namel  '_'  int2str(nj)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj+l,njmx); 
end;  end 
end 
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if  vsel==2, 

surf(X,Y,Ez(:,:,p));  shading  interp 
title(['E_z(x,y,t=',int2str(p),');  Eorward  Time  Steps: 
int2str(P)],'EontSize',  14) 

xlabel('x-axis','EontSize',  1 8);  ylabel('y-axisVEontSize',  1 8) 

axis  equal;  axis(v);  caxis(C);  colorbar 

view(0,90);  hold  on 

plot3(XBC,YBC,ZBC,’ -k’) ; 

if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,’or’);  end 

hold  off 

if  ~isempty(namel), 

if  p  ==  pj2(nj),  name=[namel  int2str(nj)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj+l,njmx); 
end;  end 
end 

if  vsel==3,  surfl(X,Y,Ez(:,:,p)) 

title(['E_z(x,y,t=',int2str(p),');  Eorward  Time  Steps:  ',int2str(P)],'PontSize',14) 
xlabel('x-axisVPontSize',  1 8);  ylabel('y-axisVPontSize',  1 8) 
axis(v);  colormap(copper); 
if  ~isempty(namel), 

if  p  ==  pjl(nj),  name=[namel  int2str(nj)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj+l,njmx); 
end;  end 
end 

Prame(p)=getframe(gcf) ; 
end 

if  vsel==l,  title([Einal  Time  Step  \pm  Log_{  10}\{IE_z(x,y,t=',... 
int2str(P),’)l\}’],’PontSize’,14) 

else 

title([Einal  Time  Step  E_z(x,y,t=',int2str(P),')'],'PontSize',14) 
end 

hcpy=input('Print  a  Hard  Copy  ?  (Y/N):  Vs'); 
if  ~isempty(hcpy),  if  hcpy  ==  'Y'  I  hcpy  ==  'y',  print;  end;  end 
name=input('Enter  "Name"  to  Save  as  Name.jpg  Pile  ?  (Hit  Enter  Key  to  Skip):  ' 
if  ~isempty(name),  eval(['print  ',name,'  -djpeg99']);  end 

PramePile=input('Enter  Pile  Name  to  Save  Movie  Prame  Array  as  *.mat 
(Press  Enter  to  Skip):  ','s'); 

if  ~isempty(PramePile),  eval(['save  ',PramePile,'  M  N  Prame']);  end 
clear  Prame 

%  Computation  Using  Reverse-Time  BC's  Stored  Prom  Sub-Grid  Boundary 

%  Constructing  Sub-Grid  Evolution  Arrays 
ds=ones(  1  ,max(Ms,N  s)) ; 

As=diag(ds(  1  :Ns- 1),  1  )-i-diag(ds(  1  :Ns- 1 ),- 1 ) ; 
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Bs=diag(ds(  1  :Ms- 1 ) ,  1  )+diag(ds(  1  :Ms- 1 ) ,- 1 ) ; 


Ps=input('Enter  time -step  to  begin  time -reversed  BC  data:  '); 

name2=input('Enter  Reverse  Soln  "Name"  for  Name_nn.jpg  ?  (Enter  Key  to  Skip):  Vs'); 
if  ~isempty(name2), 

dp2=input('Enter  Time  Step  Increment  Between  Stored  Erames:  '); 
pj2=l:dp2:Ps; 
nj=l;  njmx=length(pj2); 
end 

dispCSelect  Reverse-Time  Data  to  Use:  ') 

dispC  1  ==>  Pull  Grid  IC  and  Pull  Stored  BCs  (Exact  Reverse-Time  Solution)') 
dispC  2  ==>  No  ICs  and  Stored  BC  Subset  (Realistic  Measured  Boundary  Data)') 
Rdat=input('Make  Selection:  '); 

Ezs=zeros(Ns,Ms,Ps);  %  Reserving  3D  array  spaces 
if  Rdat  ==  1, 

%  Set  ICs  at  p=Ps  and  p=Ps-l 

Ezs(: ,  1  )=Ez(ys,xs,Ps);  Ezs(: , : ,2)=Ez(ys,xs,Ps- 1 ) ; 

end 


if  Rdat  ==  2, 

%  Explicitely  Enforce  Ez=0  BCs  on  Grid  Boundary  Por  Update 

Ezs(l,:,l)=EzBCl(:,Ps);Ezs(:,Ms,l)=EzBC2(:,Ps).'; 

Ezs(Ns, : ,  1  )=EzBC3 ( :  ,Ps) ;  Ezs( : ,  1 , 1  )=EzBC4(:  ,Ps) ; 
Ezs(l,:,2)=EzBCl(:,Ps-l);Ezs(:,Ms,2)=EzBC2(:,Ps-l).'; 

Ezs(Ns, :  ,2)=EzBC3 ( :  ,Ps- 1 ) ;  Ezs(: ,  1 ,2)=EzBC4(:  ,Ps- 1 ) . ' ; 
end 

for  p=3:Ps;  %  Equation  of  Evolution 

Ezs(:,:,p)=Ql*(As*Ezs(:,:,p-l)-i-Ezs(:,:,p-l)*Bs)-i-Q2*Ezs(:,:,p-l)-Ezs(:,:,p-2); 
%  Explicitely  Enforce  Ez=0  BCs  on  Grid  Boundary  Por  Update 
Ezs(l,:,p)=EzBCl(:,Ps-p-i-l);  Ezs(:,Ms,p)=EzBC2(:,Ps-p-i-l).'; 

Ezs(Ns, :  ,p)=EzBC3 ( :  ,Ps-p-i- 1 ) ;  Ezs( : ,  1  ,p)=EzBC4(:  ,Ps-p-i- 1 ) .' ; 
end 

elf  reset; 


ifvsel==l, 

EzScl=Ezs(:,:,l)/SPac;  %  scaling  so  abs(EzScl)  <=  DR 

[pi  ql]=fmd(0  <=  EzScl  &  EzScl  <  1);  %  nonlinear  remapping  for  lEzScll  <  1 
EzScl(pl-i-(ql-l)*Ns)=  1;  %  to  zero  log  plot 

[pi  ql]=fmd(-l  <  EzScl  &  EzScl  <  0); 

EzScl(pl-i-(ql-l)*Ns)=  -1; 

EzLog=sign(Ezs( : , : ,  1 )).  *log  1 0(abs(EzScl)) ; 
surf(Xs,Ys,EzLog);  shading  interp 
title(['\pm  Log_{  10}\{IE_z(x,y,t=',int2str(Ps),... 

')l\}  Initial  Condition  for  Time-Reversal'], 'PontSize', 14) 
xlabel('x-axis','PontSize',18);  ylabel('y-axis','PontSize',18) 
axis  equal;  axis(v);  caxis(C);  colorbar;  view(0,90);  hold  on 
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if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,’or’);  end 
hold  off 
end 


if  vsel==2, 

surf(Xs,Ys,Ezs(:,:,l));  shading  interp 
title(['E_z(x,y,t=',int2str(Ps),... 

')  Initial  Condition  for  Time-Reversal'], 'EontSize',  14) 
xlabel('x-axis','EontSize',  1 8);  ylabel('y-axis', 'EontSize',  1 8) 
axis  equal;  axis(v);  caxis(C);  colorbar;  view(0,90);  hold  on 
if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,'or');  end 
hold  off 
end 

if  vsel==3,  surfl(Ys,Xs,Ezs(:,:,l)); 

title(['E_z(x,y,t=',int2str(Ps),')  Initial  Condition  for  Time-Reversal'],'PontSize',14) 
xlabel('x-axis','PontSize',  1 8);  ylabel('y-axis','PontSize',  1 8) 
axis(v);  colormap(copper); 
end 

Prame=moviein(Ps,gcf) ; 

hcpy=input('Print  a  Hard  Copy  ?  (Y/N):  ','s'); 
if  ~isempty(hcpy),  if  hcpy  ==  'Y'  I  hcpy  ==  'y',  print;  end;  end 
name=input('Enter  "Name"  to  Save  as  Name.jpg  Pile  ?  (Hit  Enter  Key  to  Skip):  ','s') 
if  ~isempty(name),  eval(['print  ',name,'  -djpeg99']);  end 

for  p=l  :Ps;  %  Recording  animation  frames 
if  vsel==l, 

EzScl=Ezs(:,:,p)/SPac;  %  scaling  so  abs(EzScl)  <=  DR 

[pi  ql]=fmd(0  <=  EzScl  &  EzScl  <  1);  %  nonlinear  remapping  for  lEzScll  <  1 
EzScl(pl-i-(ql-l)*Ns)=  1;  %  to  zero  log  plot 

[pi  ql]=fmd(-l  <  EzScl  &  EzScl  <  0); 

EzScl(pl-i-(ql-l)*Ns)=  -1; 

EzLog=sign(Ezs( : , :  ,p)).  *log  1 0(abs(EzScl)) ; 
surf(Xs,Ys,EzLog);  shading  interp 

title(['\pm  Log_{  10}\{IE_z(x,y,t=',int2str(Ps-p-i-l),')l\};  Reverse  Time  Steps: 
int2str(Ps)] , 'EontSize'  ,14) 

xlabel('x-axis','PontSize',  1 8);  ylabel('y-axis','PontSize',  1 8) 
axis  equal;  axis(v);  caxis(C);  colorbar;  view(0,90);  hold  on 
if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,'or');  end 
hold  off 

if  ~isempty(name2), 

if  p  ==  pj2(nj),  name=[name2  '_'  int2str(nj)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj-i-l,njmx); 
end;  end 
end 

if  vsel==2, 

surf(Xs,Ys,Ezs(:,:,p));  shading  interp 
title(['E_z(x,y,t=',int2str(Ps-p-i-l),');  Reverse  Time  Steps: 
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int2stx(Ps)]  ,'FontSize' ,  1 4) 

xlabel('x-axis','FontSize',  1 8);  ylabel('y-axisVFontSize',  1 8) 
axis  equal;  axis(v);  caxis(C);  colorbar;  view(0,90);  hold  on 
if  Me  ~=  1  I  Ne  ~=  1,  plot3(XBS,YBS,ZBS,’or’);  end 
hold  off 

if  ~isenipty(nanie2), 

if  p  ==  pj2(nj),  nanie=[nanie2  int2str(nj)]; 
eval(['print  ',nanie,'  -djpeg99']);  nj=min(nj+l,njmx); 
end;  end 
end 

if  vsel==3,  surfl(Xs,Ys,Ezs(;,;,p)) 

tide(['E_z(x,y,t=',int2str(Ps-p+l),');  Reverse  Time  Steps;  ',int2str(Ps)],'PontSize',14) 
xlabel('x-axis','PontSize',  1 8);  ylabel('y-axisVPontSize',  1 8) 
v(5)=-EzMax;  v(6)=EzMax;  axis(v);  colormap(copper); 
if  ~isempty(name2), 

if  p  ==  pjl(nj),  name=[name2  int2str(nj)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj+l,njmx); 
end;  end 
end 

Prame(p)=getframe(gcf) ; 
end 

if  vsel==l,  title(['Pinal  Reversed  Time  \pm  Log_{  10}\{IE_z(x,y,t=0)l\}; 

Reverse  Time  Steps; 

int2str(Ps)]  ,'PontSize' ,  1 4) 

else, 

title(['Pinal  Reversed  Time  E_z(x,y,t=0);  Reverse  Time  Steps; 
int2str(Ps)]  ,'PontSize' ,  1 4) 
end 

hcpy=input('Print  a  Hard  Copy  ?  (Y/N);  Vs'); 
if  ~isempty(hcpy),  if  hcpy  ==  'Y'  I  hcpy  ==  'y',  print;  end;  end 
name=input('Enter  "Name"  to  Save  as  Name.jpg  Pile  ?  (Hit  Enter  Key  to  Skip);  ','s'); 
if  ~isempty(name),  eval(['print  ',name,'  -djpeg99']);  end 

PramePile=input('Enter  Pile  Name  to  Save  Movie  Prame  Array  as  *.mat 
(Press  Enter  to  Skip);  ','s'); 

if  ~isempty(PramePile),  eval(['save  ',PramePile,'  M  N  Prame']);  end 
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2.  Programs  FT_FDTD2D1  And  RT_FDTD2Dlm  For  The  Second  Example 
a.  FT_FDTD2D1  For  The  Forward  Time  Step 


clear  all 

opengl  neverselect 

%  Modified  to  use  natural  unrotated  x(cols)  y(rows)  array  definition 
Ns=input('Enter  number  of  x-points  on  sensor  sub-grid  boundary: '); 
Ms=input('Enter  number  of  y-points  in  sensor  sub-grid  boundary: '); 
N=2*Ns;M=2*Ms; 

Ez=zeros(M,N); 

x=l:N;  y=l:M;  Y=y'*ones(l,N);  X=ones(M,l)*x; 

P=input('Enter  number  of  time-steps:  '); 
namel=input('Enter  Eorward  "Name"  for  Name_nn.jpg  ? 

(Enter  Key  to  Skip):  ','s'); 
if  ~isempty(namel), 

dpl=input('Enter  Time  Step  Increment  Between  Stored  Erames:  '); 
pjl=l:dpl:P; 
nj=l;  njmx=length(pjl); 
end 

q=2;  %  q=dh/(c*dt)=2  for  half-step  algorithm 
Ql=l/(q*q);  Q2=(2-4*Q1);  cdt=l/q; 


nl=fix((N-Ns)/2)-tl;  n2=nl-tNs-l; 

ml=fix((M-Ms)/2)-i-l;  m2=ml-i-Ms-l; 

xs=nl  :n2;  ys=ml  :m2;  %  Sensor  grid  numbers 

Nb=2*(m2-ml-i-n2-nl);  %  Number  of  inner  boundary  nodes  =  2*(Ns-i-Ms-2) 

%  Defining  ordered  pairs  and  absolute  addresses  of  sub-grid  boundary  nodes 

in  (M,N)  array 

mys=zeros(Nb,l);  nxs=mys; 

mys(l  :Ns)=ml ;  mys(Ns-i- 1  :Ns-i-Ms-  l)=ml-i-l  :m2; 

mys(Ns-i-Ms:2*Ns-i-Ms-2)=m2;  mys(2*Ns-i-Ms- 1  :Nb)=m2- 1 :- 1  :ml-i-l ; 

nxs(l :Ns)=nl  :n2;  nxs(Ns+l  :Ns-i-Ms-l)=n2; 

nxs(Ns-i-Ms  :2*Ns-i-Ms-2)=n2- 1 :- 1  :nl ;  nxs(2*Ns-i-Ms- 1  :Nb)=nl ; 

%  Absolute  array  addresses  allows  no-loop  loading 
mns=(nxs-l)*M  -i-  mys; 

%  Storage  for  saving  scattered  field  boundary  data 
Ezb=zeros(Nb,P) ; 

%  User  Supplied  Eunction  Defines  Metallic  Target  Nodes  Where  Ez=0 
[my  nx]=Air2(M,N); 

Ntgt=length(nx); 

xtgt=x(nx);  ytgt=y(my);  ztgt=ones(Ntgt,l); 

mntgt=(nx-l)*M  -i-  my;  %  Absolute  array  addresses  allows  no-loop  loading 
%  Displaying  Target  Nodes  and  Subgrid 
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v=[l  N  1  M  -1.1  l.l];cv=[-l.l  1.1]; 

%  sub-grid  outline 

YBC=[ml  m2  m2  ml  ml]’;  XBC=[nl  nl  n2  n2  nl]’;  ZBC=[1  1  1  1  1]’; 

clf  reset;  surf(X,Y,Ez);  shading  interp 

axis(v);  axis  equal;  caxis(cv);  hold  on 

plot3(xtgt,ytgt,ztgt,'.k');  hold  on 

plot3(XBC,YBC,ZBC,’-k’); 

view(0,90);  figure(l) 

xlabel('x-axis','FontSize',14);  ylabel('y-axis','FontSize',14) 

title('Target  Nodes  and  Sub-grid  -  Press  a  Key  to  Continue  ...','FontSize',14); 

hold  off 

pause 

%  Constmcting  Evolution  Arrays 
d=ones(  1  ,max(N,M)) ; 

A=sparse(Q  1  *  (diag(d(  1 ;  M  - 1 ) ,  1  )-i-diag(d(  1;M-1),-1))); 

B=sparse(Q  1  *  (diag(d(  1  :N  - 1 ) ,  1  )-i-diag(d(  1;N-1),-1))); 

%  IC's  for  -X  Propagating  Gaussian  Impulse  Plane -Wave  (8  Mar  2001) 
dispCUnit  Peak  Gaussian  Impulse  Plane  Wave  Propagates  in  -x  Direction'); 
xO=input('Enter  IC  Peak  Eocation  xO  in  Grid  Units  for  the  Gaussian  Plane  Wave;  '); 
sig=input('Enter  Standard  Deviation  of  the  Gaussian  Plane  Wave  in  Grid  Units;  '); 
E0=1;  xc=(x-x0);  sig2=2*sig*sig;  ct=0;cdt:(P-l)*cdt;  xcp=x0-ct;  %  inc  wave  center 

%  Setting  IC's  at  t=-2  and  - 1  time  steps 
Ezl=Ez;  Ez2=Ez; 


%  p=-2 

Ezex=E0*ones(M,l)*exp(-((xc-2*cdt).^2)/sig2); 
Ez  1  (mntgt)=-Ezex(mntgt) ; 


%  p=-l 

Ezex=E0*ones(M,  l)*exp(  -((xc  -cdt).^2)/sig2) ; 

Ez2(mntgt)=-Ezex(mntgt) ; 

%  Using  bi-polar  log  scaling  to  retain  Ez  polarity  with  selected  dynamic  range 

DR=input('Enter  2-D  Plot  Dynamic  Range,  e.g.  100,  1000,  etc  ;  '); 

SFac=2*E0/DR;  %  Assuming  maxlEzl=2*E0 

Cmax=logl0(DR);  cv=[-Cmax  Cmax]; 

v=[l  N  1  M  -Cmax  Cmax];  ZBC=Cmax*[l  111  1]'; 

%  Computing  Scattered  Field  and  Displaying  Total  Field 

pshow=(5:5:P);  %  every  5  time-step  indices  to  display  progress 

for  p=l  :P;  %  Time-Stepping 

%  Compute  Exact  Plane  Wave 
Ezex=E0*ones(M,  1 )  *exp(  -((xc-i-ct(p))  .^2)/ sig2) ; 
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%  Equation  of  Evolution 

Ez  =  A*Ez2  +  Ez2*B  +  Q2*(Ez2-Ezl); 

%  Enforce  PEC  EC's  on  Target  Nodes 
Ez(mntgt)=-Ezex(mntgt) ; 

%  Enforce  PEC  EC's  on  Grid  Boundary 

Ez(:,l)=0;  %  x=l  BC 

Ez(l,:)=0;  %  y=l  BC 

Ez(:,N)=0;  %  x=M  BC 

Ez(M,:)=0;  %y=NBC 

%  Time -shift  arrays 

Ezl=Ez2;  Ez2=Ez; 

%  Saving  scattered  field  at  boundary 
Ezb(:  ,p)=Ez(mns) ; 

%  Displaying  +!-  LoglO(Total  Eield) 

EzScl=(Ez-i-Ezex)/SEac;  %  scaling  so  abs(EzScl)  <=  DR 

[pi  ql]=find(0  <=  EzScl  &  EzScl  <  1);  %  nonlinear  remapping  for  lEzScll  <  1 

EzScl(pl-i-(ql-l)*M)=  1;  %  to  zero  log  plot 

[pi  ql]=fmd(-l  <  EzScl  &  EzScl  <  0); 

EzScl(pl-t(ql-l)*M)= -1; 

EzLog=sign(Ez-i-Ezex).  *logl  0(abs(EzScl)) ; 

surf(X,Y,EzLog);  shading  interp 

axis(v);  axis  equal;  caxis(cv);  colorbar;  hold  on 

%  Highlight  Target  Locations  and  Sub-Grid  Boundary 

plot3(xtgt,ytgt,ztgt,'.k');  hold  on 

plot3(XBC,YBC,ZBC,’-k’); 

view(0,90); 

xlabel('x-axis',EontSize',14);  ylabel('y-axis','EontSize',14) 
title(['\pm  Log_{10}\[IE_z^{tot}(x,y,t=',mt2str(p),')l\}'],'EontSize',14); 
hold  off 
figure(l) 

if  ~isempty(namel), 

if  p  ==  pjl(nj),  name=[namel  int2str(p)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj-i-l,njmx); 
end;  end 
pause(.Ol) 
end 

name=input('Enter  "Name"  to  Save  as  Name.jpg  Eile  ?  (Hit  Enter  Key  to  Skip):  ','s'); 
if  ~isempty(name),  eval(['print  ',name,'  -djpeg99']);  end 

name2=input('Enter  "Name"  for  Name.Mat  to  Save  [N  M  P  Ns  Ms  xO  sig  nx  my  Ezb] 
(Enter  Key  to  Skip):  ','s'); 

if  ~isempty(name2),  eval(['save  ',name2,'  N  M  P  Ns  Ms  xO  sig  nx  my  Ezb']);  end 
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b. 


RT_FDTD2Dlm  For  The  Reversed  Time  Step 


clear  all 
dir  *.mat 

name=input('Enter  "Name"  for  Name.Mat  to  Retrieve  Boundary  and  Target  Data:  ','s') 
eval(['load  ',name,'  N  M  P  Ns  Ms  xO  sig  nx  my  Ezb']); 

%  N=number  of  x-points  in  fuU  grid 
%  M=number  of  y-points  in  full  grid 
%  P=number  of  time -steps 

%  Ns=number  of  x-points  on  sensor  sub-grid  boundary 
%  Ms=number  of  y-points  in  sensor  sub-grid  boundary 
%  xO=Gaussian  Plane  Wave  Initial  Peak  Location  in  Grid  Units 
%  sig=Gaussian  Plane  Wave  Standard  Deviation  in  Grid  Units 
%  nx(l,Ntgt)=x-target  nodes  in  (M,N)  full  grid 
%  my(l,Ntgt)=y-target  nodes 
%  Ezb(Nb,P)=scattered  field  at  subgrid  boundary 

nl=fix((N-Ns)/2)-tl;  n2=nl-tNs-l; 
ml=fix((M-Ms)/2)-i-l;  m2=ml-i-Ms-l; 

xs=nl  :n2;  ys=ml  :m2;  %  Inner  grid  coordinates  in  (M,N)  grid 

Ys=ys'*ones(l,Ns);  Xs=ones(Ms,l)*xs; 

Nb=2*(m2-ml-i-n2-nl);  %  Number  of  inner  boundary  nodes  =  2*(Ns-i-Ms-2) 

%  Computing  Boundary  Node  Indices  in  Inner  Sub-Grid  System 

myb=zeros(Nb,  1 ) ;  nxb=myb; 

myb(  1  :Ns)=  1 ;  myb(Ns-i- 1  :Ns-i-Ms- 1  )=2:Ms; 

myb(Ns-tMs  :2*Ns-)-Ms-2)=Ms;  myb(2*Ns+Ms- 1  :Nb)=Ms- 1 :- 1 :2; 

nxb(  1  :Ns)=  1  :Ns;  nxb(Ns-4- 1  :Ns+Ms- 1  )=Ns; 

nxb(Ns-tMs:2*Ns-tMs-2)=Ns- 1 :- 1 : 1 ;  nxb(2*Ns-tMs- 1  :Nb)=l ; 

mnb=(nxb-l)*Ms  -i-  myb;  %  Absolute  array  addresses  allows  no-loop  loading 

%  Metallic  Target  Nodes  Where  Ez=-Ez^inc 
Ntgt=length(nx);  ztgt=ones(Ntgt,l);  nn=(l:Ntgt)';  nnl=(l:Ntgt)'; 
nxl=nx;myl=my; 

dispCExact  Target  Nodes  (node#,  nx,  my)  in  Pull  Grid:') 
disp([nn  nx  my]); 

yn=input('Change  Assumed  Target  Nodes  ?  (Y/N):  ','s'); 
ifyn=='Y'lyn  ==  'y', 
while  1, 

n=input('Enter  Node#  to  Change  (0  to  End):  '); 
if  n  <  1  I  n  >  Ntgt,  break;  end 
nxmy=input('Enter  New  [nx  my]  as  Vector:  '); 
nxl(n)=nxmy(l);  myl(n)=nxmy(2); 
dispCRevised  Target  Nodes  (node#,  nx,  my)  in  Pull  Grid:') 
disp([nnl  nxl  myl]); 
end 
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end 

%  Assumed  Target  Node  Numbers  in  Inner  Grid  Coordinates 
mys=my-ml+l;  nxs=nx-nl+l; 

%  Absolute  Array  Addresses  for  Assumed  Target  Nodes 
mnstgt=(nxs-l)*Ms  +  mys; 

%  Assumed  Target  Node  Numbers  in  Inner  Grid  Coordinates 
mysl=myl-ml+l;  nxsl=nxl-nl+l; 

%  Absolute  Array  Addresses  for  Assumed  Target  Nodes 
mnstgtl=(iixsl-l)*Ms  +  mysl; 

q=2;  %  Assuming  q=dh/(c*dt)=2  for  half-step  algorithm 
Ql=l/(q*q);  Q2=(2-4*Q1);  cdt=l/q;  ct=cdt*(0:(P-l)); 

E0=1;  xcp=xO-ct;  sig2=2*sig*sig;  %  Incident  field  peak  location 


%  Constmcting  Sub-Grid  Evolution  Arrays 
ds=ones(  1  ,max(Ns,Ms)) ; 

As=Q  1  *  (diag(ds(  1 :  Ms- 1 ),  1  )-i-diag(ds(  1  :Ms- 1 ) ,- 1 ) ) ; 
Bs=Ql*(diag(ds(l:Ns-l),l)-i-diag(ds(l:Ns-l),-l)); 

Ps=input( ['Enter  time-step  <=  ',int2str(P),'  to  initiate  time -reversed  EDTD  solution:  ']); 

%  Using  bi-polar  log  scaling  to  retain  Ez  polarity  with  selected  dynamic  range 
DR=input('Enter  2-D  Plot  Dynamic  Range,  e.g.  100,  1000,  etc  :  '); 

SPac=2*E0/DR;  %  Assuming  maxlEzl=2*E0 
Cmax=logl0(DR);  cv=[-Cmax  Cmax]; 
v=[l  N  1  M  -Cmax  Cmax];  yi=[l  N]';  zi=[Cmax  Cmax]'; 
namel=input('Enter  "Name"  for  Name_nn.jpg  ?  (Enter  Key  to  Skip):  ','s'); 
if  ~isempty(namel), 

dpl=input('Enter  Time  Step  Increment  Between  Stored  Erames:  '); 
pjl=Ps:-dpl:l; 
nj=l;  njmx=length(pjl); 
end 

%  Reserving  space  for  evolution  field  arrays 
Ez=zeros(Ms,Ns);  Ezl=Ez;  Ez2=Ez; 

Energy=zeros(Ps,l);  %  Accumulated  Energy 
%  Reverse  Time  Evolution 
for  p=Ps:-l:l; 

xi=[xcp(p)  xcp(p)]';  %  Eor  Incident  Eield  Dashed  Line 
%  Initializing  BC's 
if  p==Ps  %  p=Ps 

Ezi=E0*ones(Ms,l)*exp(-((xs-xcp(Ps)).^2)/sig2);  %  Incident  field 
if  yn  ==  'Y'  I  yn  ==  'y', 

Ez  1  (mnstgt  1  )=-Ezi(mnstgt  1 ) ;  %  T arget  data 
else 
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Ezl(mnstgt)=-Ezi(mnstgt);  %  Target  data 
end 

Energy  ( 1  )=sum(sum(Ez  1 .  *Ez  1 )) ; 
elseif  p==Ps-l%  p=Ps-l 
Ez2(mnb)=Ezb(:,Ps-l);  %  Boundary  data 
Ezi=E0*ones(Ms,l)*exp(-((xs-xcp(Ps-l)).^2)/sig2);  %  Incident  field 
if  yn  ==  'Y'  I  yn  ==  'y', 

Ez2(mnstgtl)=-Ezi(mnstgtl);  %  Target  data 
else 

Ez2(mnstgt)=-Ezi(mnstgt);  %  Target  data 
end 

Energy(2)=Energy(  1  )+sum(sum(Ez2.  *Ez2)) ; 
elf  reset 
else 

Ezi=E0*ones(Ms,l)*exp(-((xs-xcp(p)).^2)/sig2);  %  Incident  Eield 

Ez  =  As*Ez2  +  Ez2*Bs  +  Q2*(Ez2-Ezl);  %  EDTD  Evolution 

if  yn==  'Y'  I  yn==  'y', 

Ez(mnstgtl)=-Ezi(mnstgtl);  %  Target  Node  PEC  BC's 

else 

Ez(mnstgt)=-Ezi(mnstgt); 

end 

Ez(mnb)=Ezb(:,p);  %  Boundary  Data  Update 

Ezl=Ez2;  Ez2=Ez;  %  Time-shift  arrays 

Energy(Ps-p-i- 1  )=Energy(Ps-p)-i-sum(sum(Ez.  *Ez)) ; 
end 

%  Displaying  +!-  LoglO(Total  Eield) 

EzScl=Ez/SEac;  %  Scattered  Eield  (scaling  so  abs(EzScl)  <=  DR) 

%  EzScl=(Ez-tEzi)/SEac;  %  Total  Eield 

[pi  ql]=find(0  <=  EzScl  &  EzScl  <  1);  %  nonlinear  remapping  for  lEzScll  <  1 

EzScl(pl-i-(ql-l)*Ms)=  1;  %  to  zero  log  plot 

[pi  ql]=find(-l  <  EzScl  &  EzScl  <  0); 

EzScl(pl-i-(ql-l)*Ms)=  -1; 

EzLog=sign(Ez).*loglO(abs(EzScl));  %  Scattered  Eield 
%  EzLog=sign(Ez-i-Ezi).*loglO(abs(EzScl));  %  Total  Eield 
surf(Xs,Ys,EzLog);  shading  interp 
axis  equal;  axis(v);  caxis(cv);  colorbar;  view(0,90);  hold  on 
%  Highlight  Target  Locations  and  Sub-Grid  Boundary 
ifyn  ==’Y’lyn==  'y'; 

plot3(nxl,myl,ztgt,'.r');  hold  on 
end 

plot3(nx,my,ztgt,'.k');  hold  on 
plot3(xi,yi,zi,'-r'); 

%  plot3(XBC,YBC,ZBC,’-k’); 

xlabel('x-axisVEontSize',14);  ylabel('y-axis','EontSize',14) 

title(['\pm  Log_{10}\{IE_z^{scat}(x,y,t=',int2str(p),')l\}'],'EontSize',14); 

hold  off 

if  ~isempty(namel), 
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if  p  ==  pjl(nj),  name=[namel  int2str(nj)]; 
eval(['print  ',name,'  -djpeg99']);  nj=min(nj+l,njmx); 
end;  end 

%figure(l) 

pause(.Ol) 

end 

nanie=input('Enter  "Name"  to  Save  as  Name.jpg  File  ?  (Hit  Enter  Key  to  Skip):  ','s') 
if  ~isempty(name),  eval(['print  ',name,'  -djpeg99']);  end 

tp=(l:Ps)’; 

figure(2);  plot(tp, Energy) 
xlabel('Reversed  Time  StepsVFontSize',14) 
ylabel('Energy',EontSize',  14) 

title([Einal  Sub-Grid  Energy  =  ',num2str(Energy(Ps))],EontSize',14) 


3.  Target  Functions 

a.  Air2:  Two  Aircraft-like  Targets 

function  [ny,  mx]=Air2(N,M) 

%  [ny  mx]  are  column  arrays  of  target  nodes  for  2  aircraft  from  RTWE_2D5.m 

%  Define  Basic  Aircraft  Shaped  Grid  Locations  Centered  at  y=x=0 
Tgt=zeros(29,2);  %  Basic  Target  Node  Indices  (Y,X)  ordered 
Tgt(l:ll,2)=(-5:5)’; 

Tgt(12,2)=-5;  Tgt(12,l)=2; 

Tgt(13,2)=-4;  Tgt(13,l)=l; 

Tgt(14,2)=-1;  Tgt(14,l)=4; 

Tgt(15,2)=-1;  Tgt(15,l)=3; 

Tgt(16,2)=  0;  Tgt(16,l)=3; 

Tgt(17,2)=  0;  Tgt(17,l)=2; 

Tgt(18,2)=  1;  Tgt(18,l)=2; 

Tgt(19,2)=  1;  Tgt(19,l)=l; 

Tgt(20,2)=  2;  Tgt(20,l)=l; 

Tgt(21:29,2)=Tgt(12:20,2); 

Tgt(21:29,l)=-Tgt(12:20,l); 

mc=fix((M-l)/2)-i-l;  nc=fix((N-l)/2)-i-l;  %  Approx  Grid  Center 

%  Defining  Nodes  for  Two  Offset  Targets 
ny=zeros(58,l);  mx=ny; 

ny(l:29)=Tgt(:,l)  -f-  nc  -i-  7;  %  y  offset  for  tgt  #1 
mx(l:29)=Tgt(:,2)  -i-  me  -i-  2;  %  x  offset  for  tgt  #1 
ny(30:58)=Tgt(:,l)  -i-  nc  -  7;  %  y  offset  for  tgt  #2 
mx(30:58)=Tgt(:,2)  -i-  me  -  2;  %  x  offset  for  tgt  #2 
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b.  Pt2:  Two  Point-Like  Targets 

function  [my,  nx]=Pt2(M,N) 

%  [my  nx]  are  column  arrays  of  simple  2-point  target  nodes 
%  27  June  01  Mod  from  Point2.m  reversing  M  and  N  roles  with  x  and  y 

Tgt=zeros(2,2);  %  Basic  Target  Node  Indices  (Y,X)  ordered 
Tgt(l:2,2)=[-5  5]’; 

mc=fix((M-l)/2)-i-l;  nc=fix((N-l)/2)-i-l;  %  Approx  Grid  Center 

%  Defining  Nodes 
my=zeros(2,l);  nx=my; 
my(l:2)=Tgt(:,l)  -i-  me; 
nx(l:2)=Tgt(:,2)  -i-  nc; 
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